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We study the ferromagnetic phase transition in a randomly layered Heisenberg magnet using large- 
scale Monte-Carlo simulations. Our results provide numerical evidence for the infinite-randomness 
scenario recently predicted within a strong-disorder renormalization group approach. Specifically, 
we investigate the finite-size scaling behavior of the magnetic susceptibility which is characterized 
by a non-universal power-law divergence in the Griffiths phase. We also study the perpendicular 
and parallel spin- wave stiffnesses in the Griffiths phase. In agreement with the theoretical predic- 
tions, the parallel stiffness is nonzero for all temperatures T < T c . In contrast, the perpendicular 
stiffness remains zero in part of the ordered phase, giving rise to anomalous elasticity. In addition, 
we calculate the in-plane correlation length which diverges already inside the disordered phase at 
a temperature significantly higher than T c . The time autocorrelation function within model A dy- 
namics displays an ultraslow logarithmic decay at criticality and a nonuniversal power-law in the 
Griffiths phase. 

PACS numbers: 75.10.Nr, 75.40.-s, 05.70.Jk 



I. INTRODUCTION 

When weak quenched disorder is added to a system un- 
dergoing a classical continuous phase transition, generi- 
cally the critical behavior will either remain unchanged 
or it will be replaced by another critical point with dif- 
ferent exponent values. Which scenario is realized de- 
pends on whether or not the clean critical point fulfills 
the Harris criterion. 1 In contrast, zero-temperature quan- 
tum phase transitions generically display much stronger 
disorder phenomena including power-law quantum Grif- 
fiths singularities, 2 4 infinite-randomness critical points 
featuring exponential instead of power-law scaling, 5,6 and 
smeared phase transitions. 7 ' 8 A recent review of these 
phenomena can be found in Ref. 9, while Ref. 10 focuses 
on metalic systems and also discusses experiments. 

The reason for the disorder effects being stronger at 
quantum phase transitions than at classical transitions 
is that quenched disorder is perfectly correlated in the 
imaginary time direction. Imaginary time behaves as an 
additional dimension at a quantum phase transition and 
becomes infinitely extended at zero temperature. There- 
fore, the impurities and defects are effectively "infinitely 
large" in this extra dimension, which makes them much 
harder to average out than the usual finite-size defects 
and so increases their influence. 

For this reason, one should also expect strong uncon- 
ventional disorder phenomena at classical thermal phase 
transitions in systems in which the disorder is perfectly 
correlated in one or more space dimensions. Indeed, such 
behavior has been observed in the McCoy- Wu model, a 
disordered classical two-dimensional Ising model having 
perfect disorder correlations in one of the two dimensions. 
In a series of papers, McCoy and Wu 11 14 showed that 
this model exhibits an unusual phase transition featuring 
a smooth specific heat while the susceptibility is infinite 
over an entire temperature range. Fisher 5,6 achieved an 
essentially complete understanding of this phase tran- 



sition with the help of a strong-disorder renormaliza- 
tion group approach (using the equivalence between the 
McCoy- Wu model and the random transverse-field Ising 
chain). He determined that the critical point is of exotic 
infinite-randomness type and is accompanied by power- 
law Griffiths singularities. In a classical Ising model with 
perfect disorder correlations in two dimensions, the dis- 
order effects are even stronger than in the McCoy- Wu 
model: the sharp critical point is destroyed, and the tran- 
sition is smeared over a range of temperatures. 15,16 

Recently, another classical system with perfect disor- 
der correlations in two dimensions was investigated by 
means of a strong-disorder renormalization group. 17 This 
theory predicts that the randomly layered Heisenberg 
magnet features a sharp critical point (in contrast to 
the Ising case discussed above). However, it is of ex- 
otic infinite-randomness type. Somewhat surprisingly, it 
is in the same universality class as the quantum critical 
point of the random transverse-field Ising chain. 

In this paper, we present the results of Monte-Carlo 
simulations of the randomly layered Heisenberg model. 
They provide numerical evidence in support of the above 
renormalization group predictions. Our paper is orga- 
nized as follows. In Sec. II, we define our model and 
discuss its phase diagram. We also briefly summarize the 
predictions of the strong disorder renormalization group 
theory. 17 In Sec. Ill, we describe our Monte-Carlo simu- 
lations, we present the results and compare them to the 
theory. We conclude in Sec. IV. 



II. MODEL AND RENORMALIZATION GROUP 
PREDICTIONS 

We consider a ferromagnet consisting of a random se- 
quence of layers made up of two different ferromagnetic 
materials, see sketch in Fig. 1. 

Its Hamiltonian, a classical Heisenberg model on a 
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FIG. 1. (Color online) Schematic of the layered Heisenberg 
magnet: It consistes of a random sequence of layers of two 
different ferromagnetic materials. 17 

three-dimensional lattice of perpendicular size L± (in z 
direction) and in-plane size L\\ (in the x and y directions) 
is given by 

h = -J24 (S r -S r+ x+S r -S r+ y)-^ jj- S r -S r+i . (1) 

r r 

Here, S r is a three-component unit vector on lattice site 
r, and x, y, and z are the unit vectors in the coordinate 
directions. The interactions within the layers, jj, and 
between the layers, J^, are both positive and indepen- 
dent random functions of the perpendicular coordinate 
z. 

In the following, we take all to be identical, = 
J x , while the jj are drawn from a binary probability 
distribution 

P(J") = (l-p)5(jW - J U )+PS(J\\ - Ji) (2) 

with J u > J\. Here, p is the concentration of the "weak" 
layers while 1 — p is the concentration of the "strong" 
layers. 

The qualitative behavior of the model (1) is easily ex- 
plained (see Fig. 2). At sufficiently high temperatures, 
the model is in a conventional paramagnetic (strongly 
disordered) phase. Below a temperature T u (the transi- 
tion temperature of a hypothetical system having jj = 
J u for all z) but above the actual critical temperature 
T c , rare thick slabs of strong layers develop local or- 
der while the bulk system is still nonmagnetic. This is 
the paramagnetic (weakly disordered) Griffiths phase (or 
Griffiths region). In the ferromagnetic (weakly ordered) 
Griffiths phase, located between T c and a temperature 
T\ (the transition temperature of a hypothetical system 
having jj = J\ for all z\ bulk magnetism coexists with 
rare nonmagnetic slabs. Finally, below TJ, all slabs are 
locally ferromagnetic and the system is in a conventional 
ferromagnetic (strongly ordered) phase. 

In Ref. 17, the behavior in both Griffiths phases and 
at criticality has been derived within a strong-disorder 
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FIG. 2. (Color online) Schematic phase diagram of the ran- 
domly layered Heisenberg magnet (1). SD and SO denote the 
conventional strongly disordered and strongly ordered phases, 
respectively. WD and WO are the weakly disordered and or- 
dered Griffiths phases. T c is the critical temperature while T u 
and Ti mark the boundaries of the Griffiths phase . 

renormalization group calculation. Here, we simply mo- 
tivate and summarize the results. The probability of find- 
ing a slab of Lrr consecutive strong layers is given by 
simple combinatorics; it reads w(Lrr) ~ (1 — p) LRR = 
e -pL RR with p — _ l n (i — p). Each such slab is equiva- 
lent to a two-dimensional Heisenberg model with an ef- 
fective interaction LrrJ u . Because the two-dimensional 
Heisenberg model is exactly at its lower critical dimen- 
sion, the renormalized distance from criticality, e, of 
such a slab decreases exponentially with its thickness, 
€(Lrr) ~ e -bL RR 9,18 Combining the two exponentials 
gives a power-law probability density of locally ordered 
slabs, 

p(e) ~ e^ 6 - 1 = e 1 ^- 1 (3) 

where the second equality defines the conventionally used 
dynamical exponent, z. It increases with decreasing tem- 
perature throughout the Griffiths phase and diverges as 
z ~ 1/\T — T c \ at the actual critical point. 

Many important observables follow from appropriate 
integrals of the density of states (3). The susceptibility 
can be estimated by x ~ / dep(e)/e. In an infinite sys- 
tem, the lower bound of the integral is 0; therefore, the 
susceptibility diverges in the entire temperature region 
where z > 1. A finite system size L\\ in the in-plane di- 
rections introduces a nonzero lower bound e m i n ~ L^ 2 . 
Thus, for z > 1, the susceptibility in the weakly disor- 
dered Griffiths phase diverges as 

X {L^)~L 2 f 2 > z (4) 

and in the weakly ordered Griffiths phase, it diverges as 

X (L||)~lJ +2/ ^. (5) 

The strong-disorder renormalization group 17 confirms 
these simple estimates and gives x ~ L 2 [\n (Ly/a)] 2 ^ -1 /^ 

at criticality where (j) = (1 + \/5)/2 and -0 = 1/2 are crit- 
ical exponents of the infinite randomness critical point. 

The spin-wave stiffness p s is defined by the work 
needed to twist the spins of two opposite boundaries by 
a relative angle 0. Specifically, in the limit of small and 
large system size, the free-energy density / depends on 
as 
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Because the randomly layered Heisenberg model is 
anisotropic, we need to distinguish the parallel spin- wave 
stiffness pi from the perpendicular spin- wave stiffness pjr . 
To calculate the parallel spin- wave stiffness, we apply 
boundary conditions at x = and x = L\\ and set L = Ly 
in Eq. (6) whereas the boundary conditions are applied 
at z = and z = L± to calculate the perpendicular spin- 
wave stiffness with L = L± in Eq. (6). 

Let us first discuss the parallel stiffness. In this case, 
the free energy difference f{0) — f(0) is simply the sum 
over all layers participating in the long-range order (each 

having the same twisted boundary conditions). Thus, pi 
is nonzero everywhere in the ordered phase. The strong- 
disorder renormalization group approach 17 predicts 

pl~m~\T-T c f (T<T C ) (7) 

where /3 = (3 — y/5)/2 is the order parameter expo- 
nent of the infinite randomness critical point. The 
parallel stiffness behaves like the total magnetization 
m = | (S r )|/(Lj_Ly), because both renormalize ad- 
ditively under the strong-disorder renormalization-group 
theory. 17 

If the twist is applied between the bottom (z = 0) 
and the top (z = L±) layers, the local twists between 
consecutive layers will vary from layer to layer. Min- 
imizing f{0) — f(0) leads to pjr ~ (1/J^y) -1 where 
J^jj are the effective couplings between the rare re- 
gions. Within the strong-disorder renormalization group 
approach, the distribution of the follows a power 

law p(Jeff) ~ ( J e"//) 1/2; " 1 - ThuS > Pa = in P art of 

the ordered Griffiths phase. It only becomes nonzero 
once z falls below 1 at a temperature T s < T c . Be- 
tween T c and T s , the system displays anomalous elas- 
ticity. Here, the free energy due to the twist scales with 
f{0) — f(0) ~ LJ 1 ~ Z . Thus, the perpendicular stiffness 
formally vanishes as pjr ~ L 1 ^ z with increasing L±. 

To study the dynamical critical behavior, a phe- 
nomenological dynamics is added to the randomly lay- 
ered Heisenberg model. The simplest case is a purely 
relaxational dynamics corresponding to model A in the 
classification of Hohenberg and Halperin. 19 

The dynamic behavior can be characterized by the av- 
erage time autocorrelation function 

C{t) = J d M S rW S r(0)>, (8) 

where S r (t) is the value of the spin at position r and time 
t. 

The behavior of C(t) in the weakly disordered Griffiths 
phase can be easily estimated. The correlation time of a 
single locally ordered slab is proportional to 1/e. 17 Sum- 
ming over all slabs using the density of states (3) then 
gives 

C(t) ~ / dep(e)e- et ^t~^ z . (9) 



The strong disorder renormalization group 
calculation 17 confirms this estimate. Moreover, at 
criticality, when z — )• oo, it gives an even slower 
logarithmic behavior 

c(t) ~ [Mt/to)}*- 1 ^. (io) 

where to is a microscopic length scale. 

III. MONTE-CARLO SIMULATIONS 

A. Overview 

In this section we report results of Monte-Carlo simu- 
lations of the randomly layered Heisenberg magnet. Be- 
cause the phase transition in this system is dominated 
by the rare regions, sufficiently large system sizes are re- 
quired in order to get reliable results. We have simulated 
system sizes ranging from L± = 90 to 800 and Lu = 10 
to 400. We have chosen J u = 1 and J t = 0.25 in Eq. (2). 
All the simulations have been performed for disorder con- 
centrations p = 0.8. With these parameter choices, the 
Griffiths region ranges from T\ w 0.63 to T u w 1.443. For 
optimal performance, we have used large numbers of dis- 
order realizations, ranging from 100 to 7200, depending 
on the system size. While studying the thermodynam- 
ics, we have used the efficient Wolff cluster algorithm 20 
to eliminate critical slowing down. We have equilibrated 
every run by 100 Monte- Carlo sweeps, and we have used 
another 100 sweeps for measurements. To investigate the 
critical dynamics, we have equilibrated the system using 
the Wolff algorithm but then propagated the system in 
time by means of the Metropolis algorithm 21 which im- 
plements model A dynamics. 

B. Thermodynamics 

To test the finite-size behavior (4, 5) of the susceptibil- 
ity, one needs to consider samples having sizes L±_ ^> L\\ 
such that L± is effectively infinite. We have used sys- 
tem sizes L± = 800 and Ln =10 to 90. Figure 3 shows 
the susceptibility x as a function of L\\ for several tem- 
peratures in the Griffiths region between T\ = 0.63 and 
T u ^ 1.443. In agreement with the theoretical predic- 
tions (4) and (5), x follows a nonuniversal power law 
in L|| with a temperature-dependent exponent. Simu- 
lations for many more temperature values, in the range 
T w 0.76 — 1.2, yield analogous results. 

The values of the exponent z extracted from fits to 
(4, 5) are shown in Fig. 4 for the paramagnetic and 
ferromagnetic sides of the Griffiths region, z can be fitted 
to the predicted power law z ~ 1/|T — T c |, as discussed 
after (3), giving the estimate T c « 0.933. 

For a deeper understanding of the thermodynamic crit- 
ical phenomena of the layered Heisenberg model, we 
have also studied the behavior of the in-plane correla- 
tion lengths in Griffiths phase. Figure 5 shows the scaled 
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FIG. 3. (Color online) Susceptibility % as a function of in- 
plane system size Ly for several temperatures in the Griffiths 
region. The perpendicular size is L± = 800; the data are 
averages over 300 disorder configurations. The solid lines are 
fits to the power laws (4, 5). 



FIG. 5. (Color online) Scaled in-plane correlation length 
£||/L|I as a function of temperature T for several in-plane sys- 
tem sizes L|| in the Griffiths region. The perpendicular size 
is L± = 800; the data are averaged over 300 disorder config- 
urations. 
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FIG. 4. (Color online) Griffiths dynamical exponent z vs 
temperature. The data are extracted from the perpendicular 
stiffness data in Fig. 6b, the susceptibility data in Fig. 3, the 
parallel correlation length data in Fig. 5 and the autocorrela- 
tion function data in Fig. 7. The solid lines are a power-law 
fit of z (extracted from Fig. 3) to (4) and (5). 



correlation length £||/^|| as a function of temperature 
for different values of Ln. Surprisingly, the curves cross 
at a temperature, T & 1.17, significantly higher than 
T c « 0.93. This implies that the average in-plane corre- 
lation length diverges in part of the disordered phase. 

To understand this behavior, we estimate the rare re- 
gion contribution to the averaged in-plane correlation 
length. It can be calculated by integrating over the den- 
sity of states (3) as 
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(ii) 



where ~ l/ e * s ^ ne dependence of the in-plane cor- 
relation length of a single region 17,22 on the renormalized 
distance e from criticality. Note that we average £ 2 in- 
stead of £|| because that is what numerically happens in 
the second moment method which defines £ 2 via 



E r C(r)r 2 



(12) 



with C(r) being the spatial correlation function. The 
integral in (11) diverges for z > 1 and converges for z < 1. 
The in-plane correlation length therefore diverges already 
in the disordered Griffiths phase at the temperature at 
which the Griffiths dynamical exponent is z = 1. From 
Fig. 5 we estimate this temperature to be T « 1.17. As 
can be seen in Fig. 4, this value is in good agreement 
with the result extracted from the finite size behavior of 
X- 

We now turn to the spin-wave stiffness. Calculating 
the stiffness by actually carrying out simulations with 
twisted boundary conditions is not very efficient. How- 
ever, the stiffness can be rewritten in terms of expecta- 
tion values calculated in a conventional run with periodic 
boundary conditions. The resulting formula which is a 
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generalization of that used by Caffarel et al 23 reads 



— ( y ] Jr,r' [S r • S r / 
W'> 



(S P .a)(S P /-a)] {z-z') 



f\2 





[(S r xS r 0-a] (z 



(13) 

Here, a can be any unit vector perpendicular to the 
total magnetization m. For pi, (z — z') has to be replaced 
by (x — x'). This formula is derived in appendix A. 

Figure 6a shows the results for the perpendicular and 
parallel stiffnesses of our randomly layered Heisenberg 
model. We have used a system of size L± = 100 and 
L|| = 400. The figure shows that the two stiffness indeed 

behave very differently. The parallel stiffness pi vanishes 
at T « 0.9 — 0.95 in good agreement with our earlier 
estimate of T c « 0.93. In contrast, the perpendicular 
stiffness vanishes at a much lower temperature T w 0.7. 
Thus, in the range between T ^ 0.7 and T c , the system 
displays anomalous elasticity, as predicted. (Note: The 
slight rounding of both pi and pjr can be attributed to 
finite-size effects.) 

The results of the perpendicular spin- wave stiffness pjr 
are analyzed in more detail in Fig. 6b for perpendicular 
sizes L± = 15 — 40. We have used a parallel size Ln = 400 
and a temperature range T = 0.65 — 0.85 where the data 
are averaged over 1000 disorder configurations. The plot 
shows a non-universal power-law dependence of pj- on 
Lj_ which agrees with the prediction 
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The dynamical exponents z extracted from fits of pjr to 
(14) are also shown in Fig. 4. While they roughly agree 
with the values extracted from x, the agreement is not 
very good. We believe this is due to the rather small L± 
values used. 



C. Critical dynamics 

To investigate the behavior of the autocorrelation func- 
tion C(t) in the weakly disordered Griffiths phase, we 
have used system sizes L± = 400 and Ln = 100 and tem- 
peratures from T = 1.25 to 1.35. From figure 7, one can 
see that the long-time behavior of C(t) in the Griffiths 
phase follows a non-universal power law which is in agree- 
ment with the prediction (9). Fits of the data to (9) can 
be used to obtain yet another estimate of the dynamical 
exponent z. The resulting values are shown in Fig. 4, 
they are in good agreement with those extracted from \. 

Figure 8 shows the behavior of C(t) near critical- 
ity plotted such that the expected logarithmic time- 
dependence (10) gives a straight line. We have used sys- 
tem sizes Lj_ = 400 and Ln = 230 and temperatures 
from T = 0.86 to 0.91. We find that C(t) indeed follows 
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FIG. 6. (Color online) a: Perpendicular and parallel spin- 
wave stiffnesses (pjr and pi , respectively) as functions of tem- 
perature T for system with sizes L± = 100 and L\\ = 400. 
The data are averaged over 100 disorder configurations, b: 
Perpendicular spin- wave stiffness as a function of L± for tem- 
peratures in the weakly ordered Griffiths phase and Ly = 400. 
The data are averaged over 1000 disorder configurations. The 
solid lines are fits to (14). 



the prediction at an estimated T c « 0.895. This esti- 
mate agrees reasonably well with that stemming from 
the finite-size behavior of x- We attribute the remaining 
difference to the finite-size effects and (in case of C(t)) 
finite-time effects. 



IV. CONCLUSIONS 

To summarize, we have reported the results of large- 
scale Monte-Carlo simulations of the thermodynamics 
and dynamic behavior of a randomly layered Heisen- 
berg model. Our results provide strong numerical evi- 
dence in support of the infinite-randomness scenario pre- 
dicted within the strong-disorder renormalization group 
approach. 17 Morever, our data are compatible with the 
prediction that the randomly layered Heisenberg model 
is in the same universality class as the one-dimensional 
random transverse-field Ising model. 
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FIG. 7. (Color online) Time autocorrelation function C(t) 
for temperatures from T — 1.25 to 1.35 (within the Griffiths 
phase). The system sizes are L± — 400 and Ly = 100. The 
data are averaged over 1720 — 7200 disorder configurations. 
The solid lines are fits to the power-law prediction (9) (with 
the fit range marked). 
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FIG. 8. (Color online) Time autocorrelation function C(t) 
for temperatures from T — 0.86 to 0.91 (near criticality) . The 
system sizes are L± = 400 and Lm = 230. The data are aver- 
aged over 70 to 80 disorder configurations. The dashed line 
shows the logarithmic behavior (10) at the estimated critical 
temperature T c = 0.895. 



An important question left unanswered by the strong- 
disorder renormalization group approach 17 is whether or 
not weakly or moderately disordered systems actually 
flow to the infinite-randomness critical point. The clean 
Heisenberg critical point is unstable against weak layered 
disorder because it violates the generalized Harris crite- 
rion d r v > 2 where d r = 1 is the number of random 
dimensions. Thus, weak layered randomness initially in- 
creases under renormalization. Our numerical parameter 
choices, p = 0.8 and J u /Ji = 4 correspond to moderate 
disorder as the distribution is not particularly broad on 
a logarithmic scale. The fact that we do confirm infinte- 
randomness behavior for these parameters suggests that 
the infinite-randomness critical point may control the 
transition for any nonzero disorder strength. A numerical 
verification of this conjecture by simulating very weakly 
disordered systems would require even larger system sizes 
and is thus beyond our present computational capabili- 
ties. 

Experimental verifications of infinite-randomness crit- 
ical behavior and the accompanying power-law Griffiths 
singularities have been hard to come by, in particular in 
higher-dimensional systems. Only very recently, promis- 
ing measurements have been reported 26 ' 27 of the quan- 
tum phase transitions in CePdi-^Rh^ and Nii-^V^. 
The randomly layered Heisenberg magnet considered 
here provides an alternative realization of an infinite- 
randomness critical point. It may be more easily real- 
izable in experiment because the critical point is classi- 
cal, and samples can be produced by depositing random 
layers of two different ferromagnetic materials. 

Magnetic multilayers with systematic variation of the 
critical temperature from layer to layer have already been 
produced, 28 and our results would apply to random ver- 
sions of these structures. 
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Appendix A: Spin-wave stiffness in terms of spin 
correlation functions 



We would have liked to determine the complete set 
of critical exponents of the infinite-randomness critical 
point directly from the numerical data. To this end we 
have attempted to perform an anisotropic finite-size scal- 
ing analysis as in Refs. 24 or 25. However, within the 
accessible range of system sizes of up to about 10 7 sites, 
the corrections to the leading scaling behavior were so 
strong that we could not complete the analysis. This 
task thus remains for the future. 



Twisted boundary conditions, i.e., forcing the spins on 
one surface of the sample of size L to make an angle of 
with those on the opposite surface, lead to a change in 
the free energy density /. It can be parametrized by 

m-m = \p.(jj ■ (ai) 

which defines the spin- wave stiffness p s . 



7 



For definiteness, assume we apply a twist of 6 around 
the perpendicular axis between the top and bottom sur- 
faces of the sample. We parametrize the Heisenberg spin 



/sin($ r ) cos(0 r )^ 
S r = sin(^ r ) sin(0 r ) 
1 cos(tf r ) 



(A2) 



The boundary conditions then read <p r = at the bottom 
(z = 0) surface and </> r = at the top (z = L±) surface. 
To eliminate the twisted boundary condition, we now 
perform the variable transformation 



(A3) 



which gives new boundary conditions of ip Y = at both 
z r = and z r = L±. 

Substituting the variable transformation in the Heisen- 
berg Hamiltonian (1), we obtain 



H 



J ry\ 

<r,r'> L 



sin($ r ) sin($ r /) 



COS ^ r — lj) T f + J^( z ~ Z*)^ + COs(^ r ) COs(^ r /)| 

(A4) 

where the twist is "distributed" over the volume. Thus, 
the twist angle now appears as a parameter of the 
Hamiltonian. We can use standard methods to refor- 
mulate the second derivative of the free energy F as 

d^F_^/dH\ 2 /d 2 H\ 1 I 'fdH\ 
30* ~ T \ 30 J + \ 80* J T \ \d0 ' ' [Ab) 



where the first term on the right hand side vanishes 
due to symmetry. Evaluating the derivatives of H for 
the Hamiltonian (A4) gives the spin-wave stiffness p s = 
L 2 (d 2 f/d0 2 )l naS 



(S r .k)(S r ,-k) (z-z') 




(S r x S r >) ■ k (z - z') 



(A6) 

Here, k is the unit vector in the z-direction. The same 
equation was derived in Ref. 23 for the XY case. Equa- 
tion A6 needs to be evaluated with fixed boundary con- 
ditions at the top and bottom layeres. Applying this for- 
mula to simulations with periodic boundary conditions 
leads to incorrect results in the Heisenberg case (even 
though it works in XY case). The reason is that Eq. 
(A6) is sensitive to twist in the XY plane only. 

In the Heisenberg case this can be fixed by aligning 
the imaginary twist axis with a direction a perpendicu- 
lar to the total magnetization in each Monte- Carlo mea- 
surement. We use a = (m x k)/|m x k|. The resulting 
formula for the spin- wave stiffness can be used efficiently 
by Monte-Carlo simulations with periodic boundary con- 
ditions. It reads 

\<r,r'} / 

~k(f[^ t ^[(S r xS r 0-a](z-^)j ^. 

(AT) 

We have tested that this equation reproduces the re- 
sults obtained directly from Eq. (Al). 
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